Homework 4

QTM 220

Summary

In this problem set, we’re going to get practice using some new regression models to do the kinds of analysis we’ve been doing all semester. While we do this, we’ll pay a bit of attention to bias—a new phenomenon, because it’s not something that happens when we use the ‘all curves’ model we’ve been using.

Here are the libraries we’ll need.

library(tidyverse)
library(rpart)

And some settings to configure how our plots and code output looks.

theme_set(theme_minimal())
options(scipen = 1, digits = 3)

Bone Density

Researchers hypothesize that higher calcium intake causes greater bone density (stronger bones). To test this, a clinic runs an experiment where patients are treated with a calcium supplement or a placebo for 1 year, after which their bone density is measured. Bone Density T-Scores indicate how good someone’s bone density is—0 is typical, positive numbers are better, negative numbers are worse (a threshold score of -2.5 is used to diagnose osteoporosis). The (fake) data from this imaginary study is shown below.

To get a random sample from its patient population in which treatment is randomized, the following procedure is used.

Patients in our population are classified as old (age 65+) or young (age 64-). We’ll call this \(X\): \(X=0\) if young, \(X=1\) if old.
For each patient in our population, we flip a coin to assign them to treatment or placebo. An age-group-specific weighted coin, so that patients are assigned treatment or placebo with probability depending on their age group. We’ll call this \(W\): \(W=1\) if treatment, \(W=0\) if placebo.
We then sample with replacement from the population.¹

Here’s a table of the data.

bone.data = data.frame(age        = c(1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0),
                       supplement = c(1, 1, 0, 1, 1, 1, 0, 0, 0, 1, 0, 1),
                       bone_density = c(-2, -3, -2, 0, -1, 0, -2, 0, 1, 1, -1, 1.5))
bone.data

age	supplement	bone_density
1	1	-2.0
1	1	-3.0
1	0	-2.0
1	1	0.0
1	1	-1.0
1	1	0.0
1	0	-2.0
1	0	0.0
0	0	1.0
0	1	1.0
0	0	-1.0
0	1	1.5

And here’s a plot. We’ll show the individual points, the mean and standard deviation of bone density for each age/treatment group (as points with arms), and the proportion of patients in each group that received treatment (as bar heights). Because we’re used to seeing data in which \(y\) is non-negative, we’ll plot bone density plus 3, so that our plot looks a bit more familiar.

plotX = function(w,x) { x + 0.15*(w-0.5) } 

bone.plot = 
  ggplot(bone.data) + 
    geom_point(aes(x=plotX(supplement, age), y=bone_density+3, color=factor(supplement)),
               alpha=1, size=2, position=position_jitter(width=0.05, height=0, seed=2)) + 
    stat_summary(aes(x=plotX(supplement, age), y=bone_density+3, color=factor(supplement)),
                 geom='pointrange', fun=mean, 
                                    fun.min=\(x)mean(x)-sd(x), 
                                    fun.max=\(x)mean(x)+sd(x), alpha=.5) +    
    geom_bar(aes(x=plotX(supplement, age), y=4*after_stat(prop), fill=factor(supplement), color=factor(supplement)), 
             position='identity', alpha=.1, linewidth=.2) + 
    scale_x_continuous(breaks=c(0,1), labels=c('Young', 'Old')) + xlab('') 
     
bone.plot

Let’s start with some exercises where we estimate population summaries by taking appropriate averages of subsample means. These should be familiar from the last problem set, although there’s more of a causal flavor to their interpretation.

Exercise

Exercise 1 Using a simple difference of column means, estimate the Conditional Average Treatment Effect, \(\CATE(x) = \mu(1,x)-\mu(0,x)\), for the young patients (i.e. for \(x=0\)) and the old patients (\(x=0)\). Here \(\mu(w,x)\) is the average bone density of patients in our population assigned to treatment \(w\) and age group \(x\).

Exercise

Exercise 2 Estimate two summaries of the Conditional Average Treatment Effect. \[ \begin{aligned} \text{ATT} &= \frac{\sum_{j:w_j=1} \CATE(x_j)}{\sum_{j:w_j=1} 1} && \text{Average Treatment Effect on the Treated Subpopulation} \\ \text{ATE} &= \frac{\sum_{j=1}^m \CATE(x_j)}{\sum_{j=1}^m 1} && \text{Average Treatment Effect on the Whole Population} \end{aligned} \]

To do this, average your estimate of the CATE from the previous exercise over the treated subsample and over the whole sample respectively.

Exercise

Exercise 3 Suppose that instead of these average treatment effect estimates, you simply took the difference in means of bone density between the treated and control groups. What is this difference in means?

Exercise

Exercise 4 What does this difference in means really estimate, i.e., what is the expected value of this difference in means? Express your answer in terms of \(\mu(w,x)\) and the proportion \(p_{w,x}=m_{w,x}/m_w\) of patients — among those assigned to treatment \(w\) — that are in age group \(x\). Is this the ATT, ATE, or something else? Is it necessarily the average of the conditional average treatment effect \(\CATE(x)\) over any group of patients? Explain why or why not.

Hint. See this slide from Lecture 11 and this slide from Lecture 12.

Now let’s recall an interpretation of what we’ve been doing Lecture 14. Using the subsample mean \(\hat\mu(w,x)\) to estimate the corresponding subpopulation mean \(\mu(w,x)\) is, equivalently, estimating the function \(\mu\) by the the function \(\hat\mu\) that minimizes the sum of squared residuals among all functions outright, i.e., \[ \hat \mu = \argmin_{m \in \model} \sum_{i=1}^n \qty{ Y_i - m(W_i,X_i) }^2 \qfor \model = \qty{ \text{ all functions } \ \ m(w,x) } \]

And we can, if we prefer, use the least squares estimator within a smaller set of functions, e.g., \[ \begin{aligned} \model &= \qty{ \text{ all functions of \emph{coarsened covariates} } \ \ m(w,x)=f(c(w,x)) } && \text{ is a piecewise-constant model } \\ \model &= \qty{ m(w,x) = m_0(w) + m_1(x) \ \ : \ \ m_0, m_1 \text{ are functions of one variable } } && \text{ is the additive model } \\ \model &= \qty{ m(w,x) = a(w) + b x \ \ : \ \ a \text{ is a function of one variable and } \ \ b \\ \text{ is a number } } && \text{ is the parallel lines model }\\ \model &= \qty{ m(w,x) = a(w) \ \ : \ \ a \text{ is a function of one variable } \ } && \text{ is the horizontal lines model }\\ \end{aligned} \]

Here’s code to find and plot you the least squares estimator \(\hat\mu\) within each model. To choose the coarsening \(c(w,x)\) in the piecewise-constant model, we’ll use a tree as implemented in the rpart package.

Here’s some generic code for making predictions and plotting given what R calls a ‘model fit object’.

prediction.function = function(model.fit) { function(w,x) { predict(model.fit, newdata=data.frame(supplement=w, age=x)) } }
plot.predictions = function(muhat) { 
  grid = expand.grid(supplement=c(0,1), age=c(0,1))
  bone.plot +
    geom_point(aes(x=age, y=muhat(supplement, age)+3, color=factor(supplement)), data=grid) + 
    geom_line(aes(x=age,  y=muhat(supplement, age)+3, color=factor(supplement)), data=grid) +
    guides(color='none', fill='none')
}

And here’s some code to actually fit the models and plot the fits. The rpart stuff is its own thing, but for the rest of them, we’re using the R built-in function lm with an appropriate regression formula describing our model. The formula language takes a bit of getting used to and we’re not going to go get into that here. But if you’re curious, feel free to come to hours and ask about it.

## The piecewise-constant model
# Tell rpart to not to worry about having more than 1 observation per level of the coarsened covariate 
# (or, in more evocative terms, per 'bucket')
rpart.settings = rpart.control(minbucket=1)  
# Then fit and plot
piecewise.constant.model.fit = rpart(bone_density ~ supplement + age, data=bone.data, control=rpart.settings)
muhat.piecewise.constant = prediction.function(piecewise.constant.model.fit)
plot.predictions(muhat.piecewise.constant) + guides(color='none')
## The additive model
additive.model.fit = lm(bone_density ~ factor(supplement) + factor(age), data=bone.data)
muhat.additive = prediction.function(additive.model.fit)
plot.predictions(muhat.additive)
## The parallel lines model
parallel.lines.model.fit = lm(bone_density ~ factor(supplement) + age, data=bone.data)
muhat.parallel = prediction.function(parallel.lines.model.fit)
plot.predictions(muhat.parallel)
## The horizontal lines model
horizontal.lines.model.fit = lm(bone_density ~ factor(supplement), data=bone.data)
muhat.horizontal = prediction.function(horizontal.lines.model.fit)
plot.predictions(muhat.horizontal)

Note

If you look at these plots closely, you can see that the least squares estimators in the additive model and the parallel lines model are the same. In fact, when \(w\) and \(x\) are binary, these two models are the same—i.e., they’re two ways of describing the same set of functions. If you like, try to prove it.

Now let’s use a few of these models and take a moment to think about what happens.

Exercise

Exercise 5 Repeat Exercise 1 - Exercise 3 using the predictions of the additive model and the horizontal lines model instead of those of the all functions model. When you do this, you’ll get identical estimates of different summaries of the population. Explain why this happens and why it didn’t happen when we used the all functions model.

Unless you actually believe these different summaries of the population are the same, this seems like a problem. And although our sample size is small, it’s not a small-sample-size problem. Using these models, we’ll get identical estimates of these different summaries no matter what sample size we have. This means at least one of the matching summary estimates has to be biased. This might be a sign that we’re using the wrong model and, in a sense, motivates our decision to use the ‘all functions’ model that we started with.

But before we move on, I want to call your attention to something. Unless you have a bug in your code (this is a good time to check), the raw difference in means you calculated in Exercise 3 and the corresponding part of Exercise 5 should be the same for every model. Let’s think about why that’s happening.

To make sense of this phenomenon, we’re going to think about the residuals we get when we predict \(Y_i\) using two functions in our model: the least squares estimator \(\hat \mu\) and a ‘perturbed version’ \(\hat\mu + t m\) for any other function \(m\) in our model. Or, to be precise, we’re going to think about the residuals we get when we predict \(Y_i\) using functions like this for a range of values of \(t\), i.e., for functions

\[ m_t(w,x) = \hat \mu(w,x) + t m(w,x) \quad \text{ for numbers } \ \ t \ \ \text{ near zero } \]

We know that the least squares estimator \(\hat \mu\) is the curve that minimizes the sum of squared residuals among all functions in our model, so we know that—if we think of the sum of squared residuals as a function of \(t\) for curves like this—the minimum occurs at \(t=0\) where \(m_t=\hat\mu\). Let’s see what the ‘zero-derivative’ condition characterizing this minimum tells us.

\[ \begin{aligned} 0 &= \frac{d}{dt} \sum_{i=1}^n \qty{ Y_i - m_t(W_i,X_i) }^2 \Big|_{t=0} && \text{ is the zero-derivative condition} \\ &= \sum_{i=1}^n 2\qty{ Y_i - m_\delta(W_i,X_i) } \times -\frac{d}{dt} m_{t}(W_i,X_i) \Big|_{t=0} && \text{ using the chain rule } \\ &= -2 \sum_{i=1}^n \qty{ Y_i - \hat\mu(W_i,X_i) } m(W_i,X_i) && \text{ differentiating } \ m_{t} \ \text{ and evaluating at } \ t=0 \\ \end{aligned} \]

This tells us that, for any function \(m\) in our model, if we multiply the least squares residuals \(\hat\varepsilon_i=Y_i-\hat\mu(W_i,X_i)\) by the corresponding values \(m(W_i,X_i)\) of that function and sum them up, we get zero. This is a generalization of the equation we used to show that the sample mean was the least squares estimator here in Lecture 14, where there was no multiplication by \(m(W_i,X_i)\)—where we just summed up the residuals.

Note

In linear algebra terms, this says that the vector \(\hat\varepsilon\) of least squares residuals is orthogonal to every vector of predictions \([m(W_1,X_1) \ldots m(W_n,X_n)]\) we can get using curves in our model. If you think of not orthogonal as meaning ‘looks a bit like’, this makes intuitive sense. If the residuals looked a bit like the predictions of some function \(m\) in our model, we’d be able to improve on the fit of \(\hat\mu\) by adding a bit of that function to it: \(\hat\mu + t m\) would, for small enough \(t\), fit the data better than \(\hat\mu\) would.

Now let’s think about what this tells us when we take \(m\) to be the constant function \(m(w,x)=1\) — assuming we’re in a model which, like every one we’ve discussed so far, includes that function. It tells us that the least squares estimator gets the overall mean right. \[ 0 = \sum_{i=1}^n \qty{ Y_i - \hat\mu(W_i,X_i) } \times 1 \quad \text{ so } \quad \frac{1}{n}\sum_{i=1}^n \hat\mu(W_i,X_i) = \frac{1}{n}\sum_{i=1}^n Y_i \]

That is, it tells us that the average of our predictions \(\hat\mu(W_i,X_i)\) over the people in our sample is the same as the average of their outcomes \(Y_i\). This is true for any model that includes the constant function \(m(w,x)=1\): the all-functions model, the additive model, the parallel lines model, the horizontal lines model, and most other models you can think of. It’s not true for the ‘lines through the origin’ model \(\{ m(w,x) = b(w)x \}\), which is why people really dislike that model.

We can say something similar about the average of our predictions \(\hat\mu(W_i,X_i)\) for people with a specific treatment status \(W_i=w\). If we use any model that includes an indicator for that treatment status, including any of the aforementioned four, then the average of our predictions \(\hat\mu(W_i,X_i)\) over the people in our sample with treatment status \(W_i=w\) is the same as the average of their outcomes \(Y_i\). \[ \begin{aligned} 0 &= \sum_{i=1}^n \qty{ Y_i - \hat\mu(W_i,X_i) } \times m(W_i,X_i) \qfor m(W_i,X_i) = \begin{cases} 1 & \text{ if } W_i=w \\ 0 & \text{ if } W_i \neq w \end{cases} \\ &= \sum_{i:W_i=w} \qty{ Y_i - \hat\mu(w,X_i) } \quad \text{ so } \frac{1}{N_w}\sum_{i:W_i=w} \hat\mu(w,X_i) = \frac{1}{N_w}\sum_{i:W_i=w} Y_i \qfor N_w=\sum_{i:W_i=w} 1 \end{aligned} \]

Exercise

Exercise 6 Use this fact to explain why the difference in means you calculated in Exercise 3 and the corresponding part of Exercise 5 is the same for every model.

Tip

Recall that the difference in means we’re talking about looks like this.

\[ \frac{1}{N_w}\sum_{i:W_i=w} \hat\mu(w,X_i) - \frac{1}{N_{w'}}\sum_{i:W_i=w'} \hat\mu(w',X_i) \quad \text{for} \quad w=1 \ \text{ and } \ w'=0 \]

What do we know about each of the terms in this difference?

Gerber, Green, and Larimer Michigan Voter Experiment

Now let’s move on to a different dataset. This one is from a field experiment conducted by Gerber, Green, and Larimer in Michigan—the same one we discussed in Lecture 6.

Let’s start by loading the data. \(Y\) will be an indicator for whether the voter turned out in the 2006 election, \(X\) be the voter’s age at the time of the election, and \(W\) be an indicator for whether the voter received the treatment. \(W=1\) indicates that they received the treatment—a letter encouraging them to turn out—and \(W=0\) indicates that they did not. To make the data a bit easier to work with, I’m going to drop voters who are 96+ years old. What’s nice about this is that it ensures we’ve got voters in both treatment groups at every remaining age.

social.pressure = read.csv('social-pressure-data.csv')
keep = social.pressure$treatment %in% c(' Control', ' Neighbors') & (2006 - social.pressure$yob)< 96
turnout.data = data.frame(W=as.numeric(social.pressure[keep,'treatment'] == ' Neighbors'),
                          X=2006-social.pressure[keep, 'yob'],
                          Y = as.numeric(social.pressure[keep,'voted'] == 'Yes'))

Here’s a plot a lot like the one we made for the bone density data. We’ll show the individual points, the mean and standard deviation of turnout for each age/treatment group (as points with arms), and the proportion of patients in each group that received treatment (as bar heights). Here the raw data is a bit harder to look at for a few reasons.

There are way more points.
There are way more columns.
The outcome is binary, so all the points in each column stack on top of each other.

To make the plot more readable, we’ll make the points transparent and jitter them a bit—both horizontally and vertically — so they’re not quite on top of each other. It helps, but it’s still a bit hard to look at the raw data. The means and proportions are still easy enough to see, though. I’m also going to plot points for a random quarter of our survey recipients, rather than all of them, to help reduce the size of the images and therefore the strain they put on our web browser. This won’t have much of an effect visually.

plotX = function(w,x) { x + 0.75*(w-0.5) } 
turnout.summary.plot = 
  ggplot(turnout.data) + 
    stat_summary(aes(x=plotX(W,X), y=Y, color=factor(W)),
                 geom='pointrange', fun=mean, 
                                    fun.min=\(x)mean(x)-sd(x), 
                                    fun.max=\(x)mean(x)+sd(x), alpha=.2, size=.5, linewidth=.5) +    
    geom_bar(aes(x=plotX(W, X), y=10*after_stat(prop), fill=factor(W), color=factor(W)), 
             position='identity', alpha=.1, linewidth=.1) + 
    guides(fill='none', color='none') + xlab('') 


plot.point = sample(1:nrow(turnout.data), round(nrow(turnout.data)/4))
turnout.plot = turnout.summary.plot + 
  geom_point(aes(x=plotX(W, X), y=Y,  color=factor(W)),
            alpha=.01, size=.01, position=position_jitter(width=0.25, height=.1, seed=2), data=turnout.data[plot.point,])  

turnout.plot

Now let’s fit a few models and plot those in. We’ll actually leave the individual points out because there are so many that the plots get a bit busy, so what we’ll be plotting over are the means and proportions. We’ll use the all-functions model, additive model, the parallel lines model, and the horizontal lines model. We already know the least squares fit of the all-functions model is going to go through the column means we’ve already plotted, so we don’t really need to plot that. But to set up later code, we’ll plot it anyway.

prediction.function = function(model.fit) { function(w,x) { predict(model.fit, newdata=data.frame(W=w, X=x)) } }
plot.predictions = function(muhat) { 
  grid = expand.grid(W=c(0,1), X=min(turnout.data$X):max(turnout.data$X))
  turnout.summary.plot +
    geom_line(aes(x=X,  y=muhat(W, X), color=factor(W)), data=grid) + guides(color='none', fill='none')
}

all.functions.model.fit = lm(Y ~ factor(W)*factor(X), data=turnout.data)
muhat.all = prediction.function(all.functions.model.fit)
plot.predictions(muhat.all)
additive.model.fit = lm(Y ~ factor(W) + factor(X), data=turnout.data)
muhat.additive = prediction.function(additive.model.fit)
plot.predictions(muhat.additive)
parallel.lines.model.fit = lm(Y ~ factor(W) + X, data=turnout.data)
muhat.parallel = prediction.function(parallel.lines.model.fit)
plot.predictions(muhat.parallel)
horizontal.lines.model.fit = lm(Y ~ factor(W),     data=turnout.data)
muhat.horizontal = prediction.function(horizontal.lines.model.fit)
plot.predictions(muhat.horizontal)

You can see that these models give us pretty different predictions at each age \(x\) and treatment status \(w\). Let’s go ahead and plot an estimate of the Conditional Average Treatment Effect \(\CATE(x) = \mu(1,x)-\mu(0,x)\) for each model, too.

x = min(turnout.data$X):max(turnout.data$X)
CATE.data = data.frame(x=x, 
                       all.functions    = muhat.all(1,x) - muhat.all(0,x),
                       additive         = muhat.additive(1,x) - muhat.additive(0,x),
                       parallel.lines   = muhat.parallel(1,x) - muhat.parallel(0,x),
                       horizontal.lines = muhat.horizontal(1,x) - muhat.horizontal(0,x)) |>
              pivot_longer(-x, names_to='model', values_to='CATE')
CATE.data$model = factor(CATE.data$model, levels=c('all.functions', 'additive', 'parallel.lines', 'horizontal.lines'))
ggplot(CATE.data) + geom_line(aes(x=x, y=CATE, color=model) ) + 
                    facet_wrap(~model, ncol=4) + 
                    ylab('Estimate of CATE(x)') + 
                    guides(color=FALSE)

Now let’s do a few exercises mirroring the major ones we did for the bone density data.

Exercise

Exercise 7 Using the four estimates of \(\CATE(x)\) plotted above, calculate corresponding estimates of the Average Treatment Effect on the Treated Subpopulation and the Average Treatment Effect on the Whole Population. \[ \begin{aligned} \text{ATT} &= \frac{\sum_{j:w_j=1} \CATE(x_j)}{\sum_{j:w_j=1} 1} && \text{Average Treatment Effect on the Treated Subpopulation} \\ \text{ATE} &= \frac{\sum_{j=1}^m \CATE(x_j)}{\sum_{j=1}^m 1} && \text{Average Treatment Effect on the Whole Population} \end{aligned} \]

To do this, average your estimate of the CATE over the treated subsample and over the whole sample respectively. You should get 8 estimates in total—an estimate of the ATT and ATE for each of the four models.

Tip

If you run the code above, then say muhat.all(1,x)-muhat.all(0,x), you’ll get an estimate of \(\CATE(x)\) based on the all-functions model. And similarly for the other models.

When you see this, you’ll notice that every estimate is the roughly the same, no matter which model you use and which summary you’re estimating. When we saw matching estimates like this before, we took it as a sign that we were using the wrong model and therefore were getting bias. But we know that can’t be what’s going on here.

Exercise

Exercise 8 Explain how you know that the estimates we get using one of these models can’t be biased.

Tip

Were we worried about bias earlier in the semester? What bearing does that have on this question?

Exercise

Exercise 9 Suppose that instead of these average treatment effect estimates, you simply took the difference in turnout predictions for the treated and control groups. That is, you calculate — letting \(\hat\mu\) be the least squares estimator of \(\mu\) based on each of your four models — this raw difference in mean predictions. \[ \hat\Delta = \frac{1}{N_1}\sum_{i:W_i=1} \hat\mu(1,X_i) - \frac{1}{N_0}\sum_{i:W_i=0} \hat\mu(0,X_i). \]

Report these estimates and compare them to the estimates of the ATT and ATE you calculated in the previous exercise. Did you really have to calculate them to know how they’d compare? Explain why or why not. You may want to refer to this slide from Lecture 12 and the proportions \(P_{w,x}\) drawn into our plots as bar heights.

Now I’m going to ask a few questions that are a bit more open-ended. They’re meant to get you thinking about the implications of these pieces. This first one is pretty straightforward.

Exercise

Exercise 10 In this study, treatment was randomized with probabilities that did not depend on age, so we knew even before looking at the data that these proportions \(P_{0,x}\) and \(P_{1,x}\), the heights of the red and green bars respectively, would be roughly the same in each column. This is enough to tell you that it wasn’t going to matter which model you used to estimate the ATT and ATE. Using your work on Exercise 9 and Exercise 6, explain why.

This second one is more challenging. I’ll make it extra credit, but I encourage you to give it a shot. You’ll probably learn something from it. And, after we review the solutions, both will be in-bounds as exam material. You’ll get a sense of what that means on the practice exam.

Exercise

Exercise 11 Extra Credit. Suppose that treatment was randomized with probabilities that depends on a coarsened version of age—the probability differs for those with \(age \le 50\) and \(age > 50\) but not within those groups. Suppose I claim that you’ll get roughly the same estimates of the ATT and ATE using the all-functions model and the not-necessarily-parallel lines model \(\model=\qty{m(w,x)=a(w)+b(w)x}\). Do you agree? If so, explain why using the same reasoning you used in Exercise 10. If not, provide a counterexample by simulating an experiment with these properties and showing that the estimates differ.